/* 

@2007, Jonathan Belof
Space Research Group
Department of Chemistry
University of South Florida

*/

#include <mc.h>


/* Gauss-Legendre integrator */
/* abscissas+weights courtesy of Abramowitz & Stegun, "Handbook of Mathematical Functions", 9th Ed., p.916 */
double rotational_integrate(system_t *system, molecule_t *molecule, int type, int li, int mi, int lj, int mj) {

	int t, p;
	double theta, phi, theta_weight, phi_weight;
	double Yi_real, Yj_real, Yi_img, Yj_img;
	double potential, integrand, integral;
	/* N = 16 */
	double roots[QUANTUM_ROTATION_GRID] =	{	-0.989400934991649932596,
							-0.944575023073232576078,
							-0.865631202387831743880,
							-0.755404408355003033895,
							-0.617876244402643748447,
							-0.458016777657227386342,
							-0.281603550779258913230,
							-0.095012509837637440185,
							 0.095012509837637440185,
							 0.281603550779258913230,
							 0.458016777657227386342,
							 0.617876244402643748447,
							 0.755404408355003033895,
							 0.865631202387831743880,
							 0.944575023073232576078,
							 0.989400934991649932596	};
	double weights[QUANTUM_ROTATION_GRID] = {	 0.027152459411754094852,
							 0.062253523938647892863,
							 0.095158511682492784810, 
							 0.124628971255533872052,
							 0.149595988816576732081, 
							 0.169156519395002538189,
							 0.182603415044923588867,
							 0.189450610455068496285,
							 0.189450610455068496285,
							 0.182603415044923588867,
							 0.169156519395002538189,
							 0.149595988816576732081,
							 0.124628971255533872052,
							 0.095158511682492784810,
							 0.062253523938647892863,
							 0.027152459411754094852	};

#ifdef XXX
	/* N = 8 */
	double roots[QUANTUM_ROTATION_GRID] = {		-0.960289856497536, -0.796666477413627, -0.525532409916329, -0.183434642495650,
						 	 0.183434642495650,  0.525532409916329,  0.796666477413627,  0.960289856497536	};
	double weights[QUANTUM_ROTATION_GRID] = {	 0.101228536290376,  0.222381034453374,  0.313706645877887,  0.362683783378362,
						 	 0.362683783378362,  0.313706645877887,  0.222381034453374,  0.101228536290376,	};

	/* N = 32 */
	double roots[QUANTUM_ROTATION_GRID] =	{	-0.997263861849481563545, -0.985611511545268335400, -0.964762255587506430774, -0.934906075937739689171,
							-0.896321155766052123965, -0.849367613732569970134, -0.794483795967942406963, -0.732182118740289680387,
							-0.663044266930215200975, -0.587715757240762329041, -0.506899908932229390024, -0.421351276130635345364,
							-0.331868602282127649780, -0.239287362252137074545, -0.144471961582796493485, -0.048307665687738316235,
							 0.048307665687738316235,  0.144471961582796493485,  0.239287362252137074545,  0.331868602282127649780,
							 0.421351276130635345364,  0.506899908932229390024,  0.587715757240762329041,  0.663044266930215200975,
							 0.732182118740289680387,  0.794483795967942406963,  0.849367613732569970134,  0.896321155766052123965,
							 0.934906075937739689171,  0.964762255587506430774,  0.985611511545268335400,  0.997263861849481563545	};

	double weights[QUANTUM_ROTATION_GRID] =	{	 0.007018610009470096600,  0.016274394730905670605,  0.025392065309262059456,  0.034273862913021433103,
							 0.042835898022226680657,  0.050998059262376176196,  0.058684093478535547145,  0.065822222776361846838,
							 0.072345794108848506225,  0.078193895787070306472,  0.083311924226946755222,  0.087652093004403811143,
							 0.091173878695763884713,  0.093844399080804565639,  0.095638720079274859419,  0.096540088514727800567,
							 0.096540088514727800567,  0.095638720079274859419,  0.093844399080804565639,  0.091173878695763884713,
							 0.087652093004403811143,  0.083311924226946755222,  0.078193895787070306472,  0.072345794108848506225,
							 0.065822222776361846838,  0.058684093478535547145,  0.050998059262376176196,  0.042835898022226680657,
							 0.034273862913021433103,  0.025392065309262059456,  0.016274394730905670605,  0.007018610009470096600	};

#endif /* XXX */


	integral = 0;
	for(p = 0; p < QUANTUM_ROTATION_GRID; p++) {
		phi = M_PI*roots[p] + M_PI;
		phi_weight = weights[p];

		for(t = 0; t < QUANTUM_ROTATION_GRID; t++) {
			theta = 0.5*M_PI*roots[t] + 0.5*M_PI;
			theta_weight = weights[t];

			/* retrieve the basis set */
			Yi_real = rotational_basis(REAL, li, mi, theta, phi);
			Yj_real = rotational_basis(REAL, lj, mj, theta, phi);
			Yi_img = rotational_basis(IMAGINARY, li, mi, theta, phi);
			Yj_img = rotational_basis(IMAGINARY, lj, mj, theta, phi);

			/* form the square modulus for integration */
			if(type == REAL)
				integrand = (Yi_real*Yj_real + Yi_img*Yj_img);
			else if(type == IMAGINARY)
				integrand = (Yi_real*Yj_img - Yi_img*Yj_real);

			/* potential part of the integrand */
			potential = molecule->quantum_rotational_potential_grid[t][p];
			integrand *= potential;

			integral += phi_weight*theta_weight*integrand;

		} /* for t */

	} /* for p */

	integral *= 0.5*M_PI*M_PI;

	return(integral);

}

